rm(list = ls())
options(warn=-1)
library(rgl)
library(misc3d)
knitr::knit_hooks$set(webgl = hook_webgl)

Load plot3D functions.

source("./plot3D_func.R")

STitch3D infers 3D spatial distributions of atrial cardiomyocytes and ventricular cardiomyocytes.

celltypes <- c("Atrial.cardiomyocytes", "Ventricular.cardiomyocytes")
celltype_colors <- c("#0000cd", "#9edae5")
um <- c(-0.8227578, -0.5278972, 0.2107004, 0,
        -0.5608635, 0.8141541, -0.1502849, 0,
        -0.0922074, -0.2418221, -0.9659296, 0,
        0, 0, 0, 1) #set the initial view of the 3D plot
spot_radius <- 0.3
axis_rescale <- c(-1,1,2)
alpha_background <- 0.1

plot3D_proportions(directory = "./results_6PCW",
                   celltypes = celltypes,
                   celltype_colors = celltype_colors,
                   um = um,
                   spot_radius = spot_radius,
                   axis_rescale = axis_rescale,
                   alpha_background = alpha_background
                   )

#the heart model is downloaded from https://github.com/MickanAsp/Developmental_heart
load('./Developmental_heart-master/data/ST_heart_spot_data.RData')
load('./Developmental_heart-master/data/heart.RData')

for (i in (0:8)){
  prop <- read.table(paste0("./results_6PCW/prop_slice", i, ".csv"), sep=",", header=TRUE)
  if (i == 0){
    prop_all <- prop
  }else{
    prop_all <- rbind(prop_all, prop)
  }
}
prop_all$idx <- unlist(lapply(strsplit(prop_all$X, "-"), function(x){x[1]}))
prop_all$sam <- unlist(lapply(strsplit(prop_all$X, "-"), function(x){substr(x[2],6,7)}))
prop_all$sam <- as.integer(prop_all$sam) + 1

atlas <- heart$atlas
atlas$idx <- atlas$spot.pos
atlas$sam <- unlist(lapply(strsplit(atlas$image, "_"), function(x){x[1]}))
atlas$sam <- as.integer(atlas$sam)

spots.table <- merge(prop_all, atlas, by=c("idx","sam"))

prop <- spots.table[, celltypes]
prop$max_prop <- apply(prop, 1, max)
prop$max_celltype <- apply(prop, 1, function(x){celltypes[which.max(x)]})

alpha_threshold = 0.2

radius.of.spots.in.atlas.pixels<- (100/(2383.36/532))/3

open3d(windowRect = c(0, 0, 720, 720))

## glX 
##   3

par3d(persp)

## NULL

view3d(userMatrix = structure(c(-0.0108720660209656, 0.899227440357208, 
                                0.437346190214157, 0, 0.955604612827301, -0.119448974728584, 
                                0.269354522228241, 0, 0.2944515645504, 0.420858442783356, 
                                -0.858007192611694, 0, 0, 0, 0, 1), .Dim = c(4L, 4L)))
drawScene.rgl(organ.dwnsmp[which(names(organ.dwnsmp)%in%c('WH'))])
for (c in 1:length(celltypes)){
  idx_celltype = (as.vector(prop$max_celltype) == celltypes[c]) &
                 (prop$max_prop > alpha_threshold)
  spheres3d(598-spots.table[idx_celltype, ]$rostral.caudal, 
            532-spots.table[idx_celltype, ]$right.left, 
            spots.table[idx_celltype, ]$anterior.posterior, 
            col=celltype_colors[c], radius=radius.of.spots.in.atlas.pixels, 
            alpha=1)
}